      program comp_disp
      implicit none
      include 'mpif.h' 
      include 'disp_header.inc'
      include 'plane_lookup.inc'
      include 'displacements.inc'
      include 'part.inc'
      character*512 parameter_file_name
      character*512 particle_base_name
      character*512 outfile_base_name
      character*512 powspec_file

      integer iargc
      integer*8 nlev
      integer nmyplanes,ifirst


      integer*8 Acoeff,Bcoeff
      integer imode
      real*8 zshift,omega0,lambda0,sigma_out
      real*8 om1
      real atemp,gfac1,disp_scale_factor,vel_scale_factor
      real lingro
      integer npmax

      integer nmaxwrite
      parameter (nmaxwrite=4)  ! Maximum number of processors able to do i/o at a time.
      integer iio, nstages
      integer*8 sum_read

      real*8 t0,t1,t2,t3
      real*8 tq0,tq1,tq2,tq3,tq4

      integer :: IERR 
      integer :: MYID
      integer :: NUMPROCS
      integer :: MYIDPLUS
      integer :: MYIDMINUS
      integer :: CYCLIC_comm

      integer*8 maxlev
      parameter (maxlev=100)
      integer*8 q_icoord1,q_icoord2
      integer*8 q_nft(maxlev),q_iwhichrand(maxlev),q_kmin(maxlev),q_kmax(maxlev)
      integer*8 q_ispheremode(maxlev), q_cic_correct(maxlev),q_iquad(maxlev)
      real*8 q_lbox(maxlev), q_downweight(maxlev),q_sigma8
      real*8 q_xcen(maxlev),q_ycen(maxlev),q_zcen(maxlev),q_sz(0:maxlev)
      integer*4 q_ns1_4(maxlev),local_nx_4, local_x_start_4,total_local_size_4
      integer*8 nlevels

!     rcrain
      integer*4, dimension(1) :: indims 
      logical, dimension(1)   :: periods = .TRUE.
      logical                 :: reorder = .TRUE.
      real                    :: xx0,yy0,zz0,rr0

      call MPI_INIT(IERR)

      ! rcrain - I added this
      call MPI_COMM_RANK(MPI_COMM_WORLD, MYID, IERR)

      call MPI_COMM_SIZE(MPI_COMM_WORLD, NUMPROCS,IERR)

      write(6,*)' NUMPROCS [post comm size] = ',NUMPROCS
      call flush(6)
      indims(:) = NUMPROCS

      call MPI_CART_CREATE(MPI_COMM_WORLD,1,indims,periods,reorder,CYCLIC_comm,IERR)

      call MPI_COMM_RANK(CYCLIC_comm, MYID, IERR)

      call MPI_CART_SHIFT(CYCLIC_comm,0,1,MYIDMINUS,MYIDPLUS,IERR)


c---------------------- I/0 pattern ---------------------------
      nstages = NUMPROCS/nmaxwrite
      if (nstages*nmaxwrite.ne.NUMPROCS) nstages=nstages+1
c--------------------------------------------------------------

      if (MYID.eq.0) then !-------------------- Processor 0 only ---
 
      if (iargc().ne.1) then
       stop 'Usage: comp_disp.x parameter_file'
      endif
      call getarg(1,parameter_file_name)

      endif ! MYID==0   ---------------------  End processor 0 only ---


c---------------------------------------
c------- Read in parameter file --------
c---------------------------------------



      call MPI_BCAST(parameter_file_name,100,MPI_CHARACTER,0,CYCLIC_comm,IERR)

      call read_paramfile(parameter_file_name,maxlev,q_nft,q_iwhichrand,
     &   q_kmin,q_kmax,q_ispheremode,q_cic_correct,q_iquad,q_lbox,q_downweight,
     &   q_xcen,q_ycen,q_zcen,q_sz,
     &   q_sigma8,zshift,omega0,lambda0,q_ns1_4,nlevels,particle_base_name,
     &   outfile_base_name,powspec_file)


      if (MYID.eq.0) t0 = MPI_WTIME()

c---------------------------------------
c-- Read in particle position data -----
c---------------------------------------


      call MPI_BCAST(particle_base_name,100,MPI_CHARACTER,0,CYCLIC_comm,IERR)
      
      if (MYID.eq.0) call check_position_files(particle_base_name)

      call MPI_BARRIER(CYCLIC_comm,ierr)

      call MPI_BCAST(ibuf_pos,12,MPI_INTEGER,0,CYCLIC_comm,IERR)

      if (MYID.eq.0) print*,'Starting to read particle data in'
      call flush(6)
      do iio=0,nstages-1
      if (mod(MYID,nstages).eq.iio) then
             call read_particles_to_processor(particle_base_name,MYID,NUMPROCS,sum_read,xx0,yy0,zz0,rr0)
      endif
      call MPI_BARRIER(CYCLIC_comm,ierr)
      if (MYID.eq.0) then
            print*,'Particle data: completed ',iio,'/',nstages
            call flush(6)
      endif
      enddo  

       call particle_read_finalize(MYID,NUMPROCS,sum_read)  
      
      if (MYID.eq.0) then
            t1 = MPI_WTIME()
            print*,'Timing: ',t1-t0,' Particle data read in and sorted '
            call flush(6)
      endif


c------------------------------------------------
c-- Calculate factors for normalising the power--
c-- spectrum normalisations and for scaling -----
c-- the velocities.                         -----
c------------------------------------------------


      atemp = 1.0/(1.0+zshift)
      gfac1 = lingro(atemp,omega0,lambda0)
      disp_scale_factor = gfac1
      om1 = (omega0*(1.+zshift)**3)/
     &     ( (omega0*(1.+zshift)**3)+(1.-omega0-lambda0)*
     &       (1.+zshift)**2 + lambda0)
      vel_scale_factor = om1**0.6   ! Approximate expression

      if (MYID.eq.0) print*,
     &   'MYID, disp_scale_factor, vel_scale_factor ',MYID,disp_scale_factor,vel_scale_factor
      call flush(6) 


c------------------------------------------------
c-- Read the power spectrum in and normalise it--
c------------------------------------------------

      if (MYID.eq.0) print*,'sigma8 = ',q_sigma8
      call flush(6)

      call fit_spline(powspec_file)   
      sigma_out = q_sigma8 * disp_scale_factor
      call normalise_powspec(sigma_out)
      call output_psp



c---------------------------------------
c-- Loop over fourier grid levels  -----
c-- and compute the displacement   -----
c-- fields                         -----
c---------------------------------------


      do nlev = 1,nlevels
c-------------------------------------------------------------------------------------
      if (MYID.eq.0) t1 = MPI_WTIME()


c  Pass params to interpolation routines

          ndimension = q_nft(nlev)
          ndim       = ndimension
          xcen       = q_xcen(nlev)
          ycen       = q_ycen(nlev)
          zcen       = q_zcen(nlev)
          sz         = q_sz(nlev)
          kmin       = q_kmin(nlev)
          kmax       = q_kmax(nlev)



c   Initialise random number generator 
          call initialise_randoms(q_iwhichrand(nlev),q_ns1_4(nlev),q_nft(nlev))

          q_icoord1         =    0     !  Calc first derivative of displ potential.

          if (MYID.eq.0) print*,'Level = ',nlev
          if (MYID.eq.0) print*,'Generating x-components ...'
          call flush(6)

          q_icoord2         =    1


          call derivs(q_icoord1,q_icoord2,dispx,q_nft(nlev),q_iwhichrand(nlev),
     &     q_lbox(nlev),q_kmin(nlev),q_kmax(nlev),q_ispheremode(nlev),
     &     q_downweight(nlev),q_cic_correct(nlev),q_iquad(nlev),
     &     local_nx_4, local_x_start_4,total_local_size_4)

          if (MYID.eq.0) print*,'Level = ',nlev
          if (MYID.eq.0) print*,'Generating y-components ...'
          call flush(6)

          q_icoord2         =    2

          call derivs(q_icoord1,q_icoord2,dispy,q_nft(nlev),q_iwhichrand(nlev),
     &     q_lbox(nlev),q_kmin(nlev),q_kmax(nlev),q_ispheremode(nlev),
     &     q_downweight(nlev),q_cic_correct(nlev),q_iquad(nlev),
     &     local_nx_4, local_x_start_4,total_local_size_4)

          if (MYID.eq.0) print*,'Level = ',nlev
          if (MYID.eq.0) print*,'Generating z-components ...'
          call flush(6)

          q_icoord2         =    3
 
          call derivs(q_icoord1,q_icoord2,dispz,q_nft(nlev),q_iwhichrand(nlev),
     &     q_lbox(nlev),q_kmin(nlev),q_kmax(nlev),q_ispheremode(nlev),
     &     q_downweight(nlev),q_cic_correct(nlev),q_iquad(nlev),
     &     local_nx_4, local_x_start_4,total_local_size_4)


         Acoeff = q_nft(nlev) * 2 * (q_nft(nlev)/2+1) !Needed for index calc
         Bcoeff =       2 * (q_nft(nlev)/2+1)         ! disp from 1-d arrays
         ifirst =  local_x_start_4                    ! in interpolation routine
         zero_point = ifirst                          ! First local fourier plane
         nmyplanes  = local_nx_4                      ! Number of planes stored locally.
        


      call flush(6)

      if (MYID.eq.0) then
             t2 = MPI_WTIME()
             print*,'Timing: Level',nlev,'Disp calc  ',t2-t1
             call flush(6)
      endif

c -----Displacements calculated - particles cycled through all processors  --------------

      call MPI_BARRIER(CYCLIC_comm,ierr)

      call MPI_ALLreduce(np,npmax,1, MPI_INTEGER, MPI_MAX, MPI_COMM_WORLD,IERR)

      call initialise_particle_interpolation_levels(MYID)

      imode = 1   !  Use trilinear interpolation

      call cycle_particles_and_interpolate(MYID,NUMPROCS, MYIDPLUS, MYIDMINUS,
     & CYCLIC_comm,npmax,Acoeff,Bcoeff,nmyplanes,imode)    


      if (MYID.eq.0) then
            t3 = MPI_WTIME()
            print*,'Timing: Level',nlev,' :',t3-t2,' to shuffle particles and compute displacements '
            call flush(6)
      endif
c-------------------------------------------------------------------------------------------------
c        If second derivatives are required for more accurate interpolation
      if (q_iquad(nlev).ne.0) then

      if (MYID.eq.0) tq0 = MPI_WTIME()


          if (MYID.eq.0) print*,'Generating dxx-components ...'
          call flush(6)

          q_icoord1         =    1     !  Calc second derivative of displ potential.
          q_icoord2         =    1


          call derivs(q_icoord1,q_icoord2,dispx,q_nft(nlev),q_iwhichrand(nlev),
     &     q_lbox(nlev),q_kmin(nlev),q_kmax(nlev),q_ispheremode(nlev),
     &     q_downweight(nlev),q_cic_correct(nlev),q_iquad(nlev),
     &     local_nx_4, local_x_start_4,total_local_size_4)

          if (MYID.eq.0) print*,'Generating dyy-components ...'
          call flush(6)

          q_icoord1         =    2
          q_icoord2         =    2

          call derivs(q_icoord1,q_icoord2,dispy,q_nft(nlev),q_iwhichrand(nlev),
     &     q_lbox(nlev),q_kmin(nlev),q_kmax(nlev),q_ispheremode(nlev),
     &     q_downweight(nlev),q_cic_correct(nlev),q_iquad(nlev),
     &     local_nx_4, local_x_start_4,total_local_size_4)

          if (MYID.eq.0) print*,'Generating dzz-components ...'
          call flush(6)

          q_icoord1         =    3
          q_icoord2         =    3
 
          call derivs(q_icoord1,q_icoord2,dispz,q_nft(nlev),q_iwhichrand(nlev),
     &     q_lbox(nlev),q_kmin(nlev),q_kmax(nlev),q_ispheremode(nlev),
     &     q_downweight(nlev),q_cic_correct(nlev),q_iquad(nlev),
     &     local_nx_4, local_x_start_4,total_local_size_4)


         Acoeff = q_nft(nlev) * 2 * (q_nft(nlev)/2+1) !Needed for index calc
         Bcoeff =       2 * (q_nft(nlev)/2+1)         ! disp from 1-d arrays
         ifirst =  local_x_start_4                    ! in interpolation routine
         zero_point = ifirst                          ! First local fourier plane
         nmyplanes  = local_nx_4                      ! Number of planes stored locally.
        
      if (MYID.eq.0) then
              tq1 = MPI_WTIME()
              print*,'Timing: dxx,dyy,dzz ',real(tq1-tq0),' Level: ',nlev
      endif

      call flush(6)

      imode = 2   !  Use quadratic interpolation on dxx,dyy,dzz components

      call cycle_particles_and_interpolate(MYID,NUMPROCS, MYIDPLUS, MYIDMINUS,
     & CYCLIC_comm,npmax,Acoeff,Bcoeff,nmyplanes,imode)  

      if (MYID.eq.0) then
              tq2 = MPI_WTIME()
              print*,'Timing: cycling part: quad disp ',real(tq2-tq1),' Level: ',nlev
      endif


          if (MYID.eq.0) print*,'Generating dyz-components ...'
          call flush(6)

          q_icoord1         =    2     
          q_icoord2         =    3


          call derivs(q_icoord1,q_icoord2,dispx,q_nft(nlev),q_iwhichrand(nlev),
     &     q_lbox(nlev),q_kmin(nlev),q_kmax(nlev),q_ispheremode(nlev),
     &     q_downweight(nlev),q_cic_correct(nlev),q_iquad(nlev),
     &     local_nx_4, local_x_start_4,total_local_size_4)

          if (MYID.eq.0) print*,'Generating dzx-components ...'
          call flush(6)

          q_icoord1         =    3
          q_icoord2         =    1

          call derivs(q_icoord1,q_icoord2,dispy,q_nft(nlev),q_iwhichrand(nlev),
     &     q_lbox(nlev),q_kmin(nlev),q_kmax(nlev),q_ispheremode(nlev),
     &     q_downweight(nlev),q_cic_correct(nlev),q_iquad(nlev),
     &     local_nx_4, local_x_start_4,total_local_size_4)

          if (MYID.eq.0) print*,'Generating dxy-components ...'
          call flush(6)

          q_icoord1         =    1
          q_icoord2         =    2
 
          call derivs(q_icoord1,q_icoord2,dispz,q_nft(nlev),q_iwhichrand(nlev),
     &     q_lbox(nlev),q_kmin(nlev),q_kmax(nlev),q_ispheremode(nlev),
     &     q_downweight(nlev),q_cic_correct(nlev),q_iquad(nlev),
     &     local_nx_4, local_x_start_4,total_local_size_4)


         Acoeff = q_nft(nlev) * 2 * (q_nft(nlev)/2+1) !Needed for index calc
         Bcoeff =       2 * (q_nft(nlev)/2+1)         ! disp from 1-d arrays
         ifirst =  local_x_start_4                    ! in interpolation routine
         zero_point = ifirst                          ! First local fourier plane
         nmyplanes  = local_nx_4                      ! Number of planes stored locally.

      if (MYID.eq.0) then
              tq3 = MPI_WTIME()
              print*,'Timing: dyz,dzx,dxy ',real(tq3-tq2),' Level: ',nlev
      endif
 
      imode = 3     ! Use quadratic interpolation of dxy,dyz,dzx components
       
      call cycle_particles_and_interpolate(MYID,NUMPROCS, MYIDPLUS, MYIDMINUS,
     & CYCLIC_comm,npmax,Acoeff,Bcoeff,nmyplanes,imode)    

      if (MYID.eq.0) then
              tq4 = MPI_WTIME()
              print*,'Timing: cycling part: quad disp ',real(tq4-tq3),' Level: ',nlev
      endif

      call flush(6)



     
      endif
c-----------------------------------------------------------------------------------------------
      enddo ! loop over levels




      if (MYID.eq.0) then
            close(10)
            t3 = MPI_WTIME()
            print*,'Timing: ',t3-t0,' excluding writing to disk '
            call flush(6)
      endif



c---------------------------------------
c-- Output to disk in stages       -----
c---------------------------------------

      if (MYID.eq.0) print*, 'Starting writing data to disk'
      call flush(6)
      do iio=0,nstages-1
       if (mod(MYID,nstages).eq.iio) then
          call write_out_binary(outfile_base_name,MYID,NUMPROCS,zshift,
     &    q_lbox(1),Omega0,Lambda0,xx0,yy0,zz0,rr0)
       endif
      call MPI_BARRIER(CYCLIC_comm,ierr)
      if (MYID.eq.0) then
         print*,'Writing ...  Completed ',iio,'/',nstages
         call flush(6)
      endif
      enddo 

      if (MYID.eq.0) then
            close(10)
            t3 = MPI_WTIME()
            print*,'Timing: ',t3-t0,' total time '
            call flush(6)
      endif



      call MPI_FINALIZE(ierr)

      end program











      subroutine check_position_files(filebase)
      implicit none
      include 'part.inc'
      character*(*) filebase
      character*512 filename
      integer numpart,if,nf
      integer*8 ntotal,itot
      integer j

      call file_name_generate(filebase,0,filename)

    
      open (12,file=filename,form='unformatted',status='old')
      read (12) ibuf_pos
      close(12)

      ntotal = nlist

      write(6,*)' *** ibuf_pos  = ',ibuf_pos
      write(6,*)' *** nfile_p   = ',nfile_p
      write(6,*)' *** filebase  = ',filebase
      write(6,*)' *** filename  = ',filename
      call flush(6)

      do j=1,nfile_p-1

      call file_name_generate(filebase,j,filename)
    
      open (12,file=filename,form='unformatted',status='old')
      read (12) numpart,itot,if,nf 
      close(12)
      write(6,*)' *** itot = ',itot
      call flush(6)
      if (itot.ne.itotal) stop 'Error:  itot and itotal do not match'
      if (j.ne.if) stop 'Error: mismatch in j,if '
      if (nfile_p.ne.nf) stop 'Error: nfile_p and nf mismatch '
      ntotal = ntotal + numpart
      if ((j.lt.nfile_p-1).and.(numpart.ne.nlist)) stop 'Error in numpart/nlist '

      enddo

      if (itotal.ne.ntotal) then
         write(6,*)'itotal = ',itotal
         write(6,*)'ntotal = ',ntotal
      endif
      if (itotal.ne.ntotal) stop 'Error: with number of particles listed '

 
     

      return
      end







      subroutine file_name_generate(filebase,ifile,filename)
      implicit none
      character*(*) filebase,filename
      integer ifile
      character*10 numeral
      integer len_trim

      write (numeral,'(i10)') ifile
      filename = filebase(1:len_trim(filebase))//'.'//
     &  numeral(1+scan(numeral,' ',.true.):len_trim(numeral))
      return
      end
      


      subroutine abort(message)
      implicit none
      character*(*) message

      print*,message
      call flush(6)
      stop
      end



      subroutine read_particles_to_processor(filebase,MYID,NUMPROCS,sum_read,xx0,yy0,zz0,rr0)
      implicit none
      include 'mpif.h'
      include 'part.inc'
      character*(*) filebase
      character*512 filename
      integer MYID, NUMPROCS

      integer nsample
      parameter (nsample=1)

      integer ndiv,nrem,nskip
      integer*8 nstart,nend
      integer ndiv_f, nrem_f,ifile1,ifile2,nind1,nind2

      integer iex,i,j
      integer iwarn 
      integer*8 sum_read
     
      real xx0,yy0,zz0,rr0

      xx0 = 0.0e0
      yy0 = 0.0e0
      zz0 = 0.0e0
      rr0 = 0.0e0

      sum_read = 0

      print*,'itotal,nsample,NUMPROCS = ',itotal,nsample,NUMPROCS
      ndiv = itotal/(nsample*NUMPROCS)
      print*,'ndiv = ',ndiv
      call flush(6)

      nrem = itotal - (nsample*NUMPROCS-1)*ndiv

      ndiv_f = itotal/nfile_p
      nrem_f = itotal - (nfile_p-1)*ndiv_f

      if (MYID.ge.NUMPROCS-2) then
      print*, 'ndiv,nrem ',ndiv,nrem
      print*, 'ndiv_f,nrem_f',ndiv_f,nrem_f
      call flush(6)
      endif



      do i=0,nsample-1
        nstart = ndiv * (NUMPROCS * i + mod(MYID,NUMPROCS))

        if ((MYID.eq.NUMPROCS-1).and.(i.eq.nsample-1)) then
           iex = 1
        else
           iex = 0
        endif
 
        if (iex.eq.1) then
                nend   = nstart + nrem - 1
        else
                nend   = nstart + ndiv - 1
        endif

        ifile1 = nstart/ndiv_f
        ifile2 = min(nend/ndiv_f,nfile_p-1)

      if (MYID.ge.NUMPROCS-2) then
      print*, 'ifile1,ifile2 ',ifile1,ifile2
      call flush(6)
      endif
         
        do j=ifile1,ifile2
        nind1 = max(nstart - j*ndiv_f,0)
        if ((iex.eq.1).and.(j.eq.nfile_p-1)) then 
              nind2 = min(nend-j*ndiv_f,nrem_f-1)
        else
              nind2 = min(nend-j*ndiv_f,ndiv_f-1)
        endif

        call file_name_generate(filebase,j,filename)

        if (j.lt.nfile_p-1) then
              nskip = nlist
        else
              nskip = nrem_f
        endif

      ! rcrain -- I've added this enable reading of the (x,y,z,r) values
      ! required to be passed into the ICs for TRAPHIC
      open (10,file=filename,form='unformatted',status='unknown')
      read (10)
      read (10) xx0,yy0,zz0,rr0
      close(10)
      call flush(6)

        call open_file(filename,iwarn) 
        
        if (MYID.ge.NUMPROCS-2) then
      print*,'filename: ',filename(1:len_trim(filename))
      print*, 'nskip,nind1,nind2 ',nskip,nind1,nind2
      call flush(6)
      endif
      
    

      call read_position_info(nskip,nind1,nind2,x(sum_read),y(sum_read),
     &    z(sum_read),ms(sum_read)  )       
        call close_file()

        if(myid == 0)then
           write(6,*)' XXX1 ',x(0),y(0),z(0)
        endif
        if(myid == 0)then
           write(6,*)' XXX2 ',x(0),y(0),z(0)
        endif
        if(myid == 0)then
           write(6,*)' XXX3 ',x(0),y(0),z(0)
        endif
        if(myid == 0)then
           write(6,*)' XXX4 ',x(0),y(0),z(0)
        endif

        sum_read = sum_read + (nind2-nind1+1)
        enddo



      enddo


      return
      end


      subroutine particle_read_finalize(MYID,NUMPROCS,sum_read) 
      implicit none
      include 'mpif.h'
      include 'part.inc'

      integer MYID, NUMPROCS
      integer IERR
      integer i,j
      integer*8 sum_read,total_read
      real*8 msum    
     

      call MPI_ALLreduce(sum_read,total_read,1,MPI_INTEGER8,
     &    MPI_SUM, MPI_COMM_WORLD,IERR)


      if (total_read.ne.itotal) stop 'Error: number of particles read in is wrong'


      msum = 0.0
      do i=0,sum_read-1
       msum = msum + ms(i)
       do j=1,3
        dxx(j,i) = 0.0  ! Set initial displacements to zero. 
       enddo
      enddo

      print*,'MYID, msum, NUMPROCS,sum_read ', MYID, msum, NUMPROCS,sum_read
      call flush(6)

      msum = msum * real(NUMPROCS)

      if (abs(msum-1.0).gt.0.05) then
         print*,'WARNING: possible load inbalance MYID, msum*NUMPROCS ',msum
         call flush(6)
      endif 


      np = sum_read    ! Number of particles read in by local processor
 

      return
      end



      subroutine initialise_particle_interpolation_levels(MYID)
      implicit none
      include 'part.inc'
      include 'disp_header.inc'

      integer MYID

      integer*4 i,ninter
      integer*8 ninter8
      real*8 dx,dy,dz,lmax
      integer nhigh, nlow
      integer l
      integer lmin
      integer ndim_i4

      nhigh = 0
      nlow  = 0
      ndim_i4 = ndimension
 
      do i=0,np-1
        ninter8 = 1+int(1.95*real(kmax)*(ms(i)**0.3333333333)/sz)
        if (ninter8.gt.32767) then
           ninter = 32767
        else
           ninter = ninter8
        endif 


        call get_cell_displacements(x(i),y(i),z(i),xcen,ycen,zcen,dx,dy,dz,lmax)

        if (lmax.lt.0.75d0) then

c   Particle lies with displacement field grid.


        nhigh = nhigh+1

        l = dble(ndimension)*dx

        lmin = mod(ndim_i4 + l - ninter + 1, ndim_i4)

        if (lmin.gt.32767) stop 'Error in lmin value - too large for i*2 '

        ninterp(i) = ninter    ! Record interpolation level 


        else
        ninterp(i) = 0         !  No interpolation required
        nlow = nlow +1
        endif

      enddo



      return
      end



      subroutine cycle_particles_and_interpolate(MYID,NUMPROCS, MYIDPLUS, MYIDMINUS,
     &   CYCLIC_comm,npmax,Acoeff,Bcoeff,nmyplanes,imode)
      implicit none
      include 'mpif.h' 
      include 'part.inc'
      integer MYID,NUMPROCS, MYIDPLUS, MYIDMINUS,CYCLIC_comm
      integer npmax
      integer*8 Acoeff,Bcoeff
      integer nmyplanes,imode

      integer i
      integer ishift,stag,rtag,lsend,IERR
      integer status(MPI_STATUS_SIZE)
      integer ntimes,nmaxrem
      integer nmaxtrans
      parameter (nmaxtrans = 50000)
      integer ijk,ipos

      do i=0,NUMPROCS-1 
  
      ishift = npmax           ! Number of particles to shift

      if (imode.eq.1) then
           call interpolate_displacements_trilinear(nmyplanes,Acoeff,Bcoeff)
      else
           call interpolate_displacements_quadratic(nmyplanes,Acoeff,Bcoeff,imode)
      endif


      stag = 0
      rtag = 0



      ntimes = ishift/nmaxtrans
      nmaxrem = ishift - ntimes*nmaxtrans
      if (nmaxrem.gt.0) then
           ntimes = ntimes + 1
      else
          nmaxrem = nmaxtrans
      endif

      if (MYID.eq.0) then
            print*,'Sending particles:',i
            print*,'ishift,ntimes,nmaxtrans,nmaxrem',ishift,ntimes,nmaxtrans,nmaxrem
            call flush(6)
      endif


      call MPI_Sendrecv_replace(np,1,MPI_INTEGER,MYIDMINUS,stag,MYIDPLUS,
     &      rtag, CYCLIC_comm,STATUS, IERR) 

      ipos = 0
       do ijk=1,ntimes
         
         if (ijk.lt.ntimes) then
            lsend = nmaxtrans
         else
            lsend = nmaxrem
         endif
       
         call MPI_Sendrecv_replace(x(ipos),lsend,MPI_REAL8,MYIDMINUS,stag,MYIDPLUS,
     &      rtag, CYCLIC_comm,STATUS, IERR)
         call MPI_Sendrecv_replace(y(ipos),lsend,MPI_REAL8,MYIDMINUS,stag,MYIDPLUS,
     &      rtag, CYCLIC_comm,STATUS, IERR)
         call MPI_Sendrecv_replace(z(ipos),lsend,MPI_REAL8,MYIDMINUS,stag,MYIDPLUS,
     &      rtag, CYCLIC_comm,STATUS, IERR)
c         call MPI_Sendrecv_replace(dxx(1,ipos),3*lsend,MPI_REAL,MYIDMINUS,stag,MYIDPLUS,
c     &      rtag, CYCLIC_comm,STATUS, IERR)
         call MPI_Sendrecv_replace(idxx(3*ipos),3*lsend,MPI_INTEGER,MYIDMINUS,stag,MYIDPLUS,
     &      rtag, CYCLIC_comm,STATUS, IERR)
         call MPI_Sendrecv_replace(ms(ipos),lsend,MPI_REAL,MYIDMINUS,stag,MYIDPLUS,
     &      rtag, CYCLIC_comm,STATUS, IERR)
         call MPI_Sendrecv_replace(ninterp(ipos),lsend,MPI_INTEGER2,MYIDMINUS,stag,MYIDPLUS,
     &      rtag, CYCLIC_comm,STATUS, IERR)


          ipos = ipos + lsend

          enddo  ! loop over sends





      enddo  !loop over processor number 


      end

      
      subroutine interpolate_displacements_trilinear(nmyplanes,Acoeff,Bcoeff)
      implicit none
      include 'part.inc'
      include 'disp_header.inc'
      include 'plane_lookup.inc'
      include 'displacements.inc'
      integer nmyplanes
      integer*8 i,j,k
      integer*8 l,m,n
      integer*8 Acoeff,Bcoeff
      integer maxplanes
      parameter (maxplanes=32767)
      integer nplane_list(maxplanes),wplane(maxplanes),mplane
      real*8 dx,dy,dz,lmax
      real*8 pweight,xweight,nroot3,yweight,zweight,tweight
      integer iymin,izmin,iy,iz
      real*8 cumulo_x,cumulo_y,cumulo_z
      integer*8 ipt
      integer*8 neff
      real*8 delta_x,delta_y,delta_z
       



      do i=0,np-1  ! loop over particles

       if (ninterp(i).gt.0) then

       neff = ninterp(i) 

       call get_cell_displacements(x(i),y(i),z(i),xcen,ycen,zcen,dx,dy,dz,lmax)


        if (lmax.lt.0.75d0) then


        if (lmax.lt.0.5) then
           pweight = 1.0d0
        else
           pweight = max(0.0d0,4.0d0*(0.75d0-lmax))
        endif
        
            
c------------- Interpolate displacements to particle position dx,dy,dz

        l = dble(ndimension)*dx
        m = dble(ndimension)*dy
        n = dble(ndimension)*dz

        delta_x = 0.5 + (dble(ndimension)*dx - dble(l) - 0.5d0)/dble(neff)
        delta_y = 0.5 + (dble(ndimension)*dy - dble(m) - 0.5d0)/dble(neff)
        delta_z = 0.5 + (dble(ndimension)*dz - dble(n) - 0.5d0)/dble(neff)
       

        call which_planes(l,ninterp(i),nmyplanes,nplane_list,wplane,mplane) 


        if (mplane.gt.0) then

         nroot3 = 1.0d0/dble(ninterp(i))
        

 
         do j=1,mplane         

         cumulo_x = 0.0   ! accumulate over individual planes
         cumulo_y = 0.0
         cumulo_z = 0.0

         if (wplane(j).gt.0) then
             xweight = delta_x * nroot3
         else
             xweight = (1.-delta_x) * nroot3
         endif

             k = nplane_list(j)

             do iymin = m-ninterp(i)+1,m + ninterp(i)     

              if (iymin.gt.m) then
                 yweight = delta_y * nroot3
              else
                 yweight = (1.-delta_y)*nroot3
              endif
 
              iy = mod(iymin+ndim,ndim)
              do izmin = n-ninterp(i)+1, n + ninterp(i)    
               if (izmin.gt.n) then
                 zweight = delta_z * nroot3
               else
                 zweight = (1.-delta_z)*nroot3
               endif

                iz = mod(izmin+ndim,ndim)

                ipt = Acoeff * k + Bcoeff * iy + iz
                tweight = xweight*yweight*zweight

                cumulo_x = cumulo_x + dispx(ipt) * tweight
                cumulo_y = cumulo_y + dispy(ipt) * tweight
                cumulo_z = cumulo_z + dispz(ipt) * tweight


              enddo
             enddo

             dxx(1,i) = dxx(1,i) + pweight*sz*cumulo_x
             dxx(2,i) = dxx(2,i) + pweight*sz*cumulo_y
             dxx(3,i) = dxx(3,i) + pweight*sz*cumulo_z

         enddo  ! loop over planes

        endif

c------------- end particle interpolation ------------------------------
        endif
       endif
      enddo    ! loop over particles

      return
      end



      subroutine interpolate_displacements_quadratic(nmyplanes,Acoeff,Bcoeff,imode)
      implicit none
      include 'part.inc'
      include 'disp_header.inc'
      include 'plane_lookup.inc'
      include 'displacements.inc'
      integer nmyplanes
      integer*8 i,j,k
      integer*8 l,m,n
      integer*8 Acoeff,Bcoeff
      integer imode
      integer maxplanes
      real*8 PI
      parameter (PI=3.1415926535897d0)
      parameter (maxplanes=32767)
      integer nplane_list(maxplanes),wplane(maxplanes),mplane
      real*8 dx,dy,dz,lmax
      real*8 pweight,xweight,nroot3,yweight,zweight,tweight
      integer iymin,izmin,iy,iz
      real*8 cumulo_x,cumulo_y,cumulo_z
      integer*8 ipt
      integer*8 neff
      real*8 delta_x,delta_y,delta_z
      real*8 rdelta
      real*8 rsignx,rsigny,rsignz
       
      if ((imode.lt.2).or.(imode.gt.3)) stop 'Error in imode value!'
    
      rdelta = -0.5d0/dble(ndim)

      do i=0,np-1  ! loop over particles

       if (ninterp(i).gt.0) then

       neff = ninterp(i) 

       call get_cell_displacements(x(i),y(i),z(i),xcen,ycen,zcen,dx,dy,dz,lmax)


        if (lmax.lt.0.75d0) then


        if (lmax.lt.0.5) then
           pweight = 1.0d0
        else
           pweight = max(0.0d0,4.0d0*(0.75d0-lmax))
        endif
        
            
c------------- Interpolate displacements to particle position dx,dy,dz

        l = dble(ndimension)*dx
        m = dble(ndimension)*dy
        n = dble(ndimension)*dz

        delta_x = 0.5 + (dble(ndimension)*dx - dble(l) - 0.5d0)/dble(neff)
        delta_y = 0.5 + (dble(ndimension)*dy - dble(m) - 0.5d0)/dble(neff)
        delta_z = 0.5 + (dble(ndimension)*dz - dble(n) - 0.5d0)/dble(neff)
       

        call which_planes(l,ninterp(i),nmyplanes,nplane_list,wplane,mplane) 


        if (mplane.gt.0) then

         nroot3 = 1.0d0/dble(ninterp(i))
        

 
         do j=1,mplane         

         cumulo_x = 0.0   ! accumulate over individual planes
         cumulo_y = 0.0
         cumulo_z = 0.0

         if (wplane(j).gt.0) then
             xweight = delta_x * nroot3
             rsignx  =  1.0d0
         else
             xweight = (1.-delta_x) * nroot3
             rsignx  = -1.0d0
         endif

             k = nplane_list(j)

             do iymin = m-ninterp(i)+1,m + ninterp(i)     

              if (iymin.gt.m) then
                 yweight = delta_y * nroot3
                 rsigny  = 1.0d0
              else
                 yweight = (1.-delta_y)*nroot3
                 rsigny = -1.0d0
              endif
 
              iy = mod(iymin+ndim,ndim)
              do izmin = n-ninterp(i)+1, n + ninterp(i)    
               if (izmin.gt.n) then
                 zweight = delta_z * nroot3
                 rsignz  = 1.0d0
               else
                 zweight = (1.-delta_z)*nroot3
                 rsignz = -1.0d0
               endif

                iz = mod(izmin+ndim,ndim)

                ipt = Acoeff * k + Bcoeff * iy + iz

           if (imode.eq.2) then  !  ie  dxx,dyy,dzz component
 
                cumulo_x = cumulo_x + dispx(ipt) * rsignx*rdelta *yweight*zweight*
     &                                                        xweight*(1.0d0-xweight)
                cumulo_y = cumulo_y + dispy(ipt) * rsigny*rdelta *zweight*xweight*
     &                                                        yweight*(1.0d0-yweight)
                cumulo_z = cumulo_z + dispz(ipt) * rsignz*rdelta *xweight*yweight*
     &                                                        zweight*(1.0d0-zweight)
          else                  ! ie dxy,dyz,dzx components

 
                cumulo_x = cumulo_x + dispy(ipt) * rsignz*rdelta *xweight*yweight*
     &                                                        zweight*(1.0d0-zweight) 
     &                              + dispz(ipt) * rsigny*rdelta *xweight*zweight*
     &                                                        yweight*(1.0d0-yweight) 


                cumulo_y = cumulo_y + dispz(ipt) * rsignx*rdelta *yweight*zweight*
     &                                                        xweight*(1.0d0-xweight) 
     &                              + dispx(ipt) * rsignz*rdelta *yweight*xweight*
     &                                                        zweight*(1.0d0-zweight) 

                cumulo_z = cumulo_z + dispx(ipt) * rsigny*rdelta *zweight*xweight*
     &                                                        yweight*(1.0d0-yweight) 
     &                              + dispy(ipt) * rsignx*rdelta *zweight*yweight*
     &                                                        xweight*(1.0d0-xweight) 


          endif

              enddo
             enddo

             dxx(1,i) = dxx(1,i) + pweight*sz*cumulo_x
             dxx(2,i) = dxx(2,i) + pweight*sz*cumulo_y
             dxx(3,i) = dxx(3,i) + pweight*sz*cumulo_z

         enddo  ! loop over planes

        endif

c------------- end particle interpolation ------------------------------
        endif
       endif
      enddo    ! loop over particles

      return
      end
      

      subroutine which_planes(l,ninterp,nmyplanes,nplane_list,wplane,mplane)
      implicit none
      include 'plane_lookup.inc'
      integer*2 ninterp
      integer*8 l
      integer nmyplanes,nplane_list(*),wplane(*),mplane
      integer j,icount
      integer lref


      lref = mod(int(l-zero_point+ndim),ndim)

c  Check if particle does not require any plane currently on the processor 


      if ((lref-ninterp.gt.nmyplanes).and.(ndim-lref.gt.ninterp)) then
            mplane = 0
            return
      endif



      icount = 0

 
      do j= l-ninterp+1,l+ninterp
        lref = mod(j-zero_point+ndim,ndim)
        if ((lref.ge.0).and.(lref.lt.nmyplanes)) then
           icount = icount + 1
           nplane_list(icount) = lref
           wplane(icount) = j-l
        endif
      enddo

      mplane = icount
      return
      end


      subroutine get_cell_displacements(xp,yp,zp,xc,yc,zc,dx,dy,dz,lmax)
      implicit none
      include 'disp_header.inc'

      real*8 xp,yp,zp,xc,yc,zc,dx,dy,dz,lmax
        
      dx = xp-xc
      dy = yp-yc
      dz = zp-zc

      if (dx.ge. 0.5d0)  dx  = 1.0d0-dx
      if (dx.lt.-0.5d0)  dx  = 1.0d0+dx       

      if (dy.ge. 0.5d0)  dy  = 1.0d0-dy
      if (dy.lt.-0.5d0)  dy  = 1.0d0+dy       

      if (dz.ge. 0.5d0)  dz  = 1.0d0-dz
      if (dz.lt.-0.5d0)  dz  = 1.0d0+dz

      dx = dx/sz 
      dy = dy/sz 
      dz = dz/sz

      lmax = max(abs(dx),abs(dy),abs(dz))  ! Distance from cube centre in units of sidelengths

      dx = dx + 0.5d0
      dy = dy + 0.5d0 
      dz = dz + 0.5d0

      dx = dx - floor(dx)  !  Give position in primary cube. Use lmax to decide
      dy = dy - floor(dy)  !  whether particle is inside the cube or not.
      dz = dz - floor(dz)  !  lmax<0.5 implies particle is within cube.

      return
      end



      subroutine write_out_binary(outfile_base,MYID,NUMPROCS,zshift,
     &   BoxSize,Omega0,Lambda0,xx0,yy0,zz0,rr0)
      implicit none
      include 'part.inc'
      include 'disp_header.inc'
      character*(*) outfile_base
      integer MYID, NUMPROCS
      real*8 zshift,BoxSize,Omega0,Lambda0
      integer*8 first_index 
      integer i
      character*512 outfile
      real xx0,yy0,zz0,rr0

c   Gadget header info:


       integer*4 unused(9) 

    

      if (MYID.ne.NUMPROCS-1) then
         first_index = MYID*np+1
      else
         first_index =itotal-np+1
      endif



      call file_name_generate(outfile_base,MYID,outfile)

      write(6,*)trim(adjustl(outfile))
      write(6,*)xx0,yy0,zz0,rr0

      open (15,file=outfile,form='unformatted',status='unknown')
      write(15) np, BoxSize, NUMPROCS, MYID, first_index, zshift,Omega0,lambda0,
     &    itotal,unused         
      write(15)xx0,yy0,zz0,rr0

      write (15) (x(i),i=0,np-1)
      write (15) (y(i),i=0,np-1)
      write (15) (z(i),i=0,np-1)
      write (15) (dxx(1,i),i=0,np-1)
      write (15) (dxx(2,i),i=0,np-1)
      write (15) (dxx(3,i),i=0,np-1)
      write (15) (ms(i),i=0,np-1)
 
      close(15)

      return
      end







 
